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It is widely accepted that the Poisson-Boltzmann (PB) theory provides a valid description for 
charged surfaces in the so-called weak coupling limit. Here, we show that the image charge repulsion 
creates a depletion boundary layer that cannot be captured by a regular perturbation approach. 
The correct weak-coupling theory must include the self-energy of the ion due to the image charge 
interaction. The image force qualitatively alters the double layer structure and properties, and 
gives rise to many non-PB effects, such as nonmonotonic dependence of the surface energy on 
concentration and charge inversion. In the presence of dielectric discontinuity, there is no limiting 
condition for which the PB theory is valid. 


I. INTRODUCTION 

The electric double layer resulting from a charged sur¬ 
face in an aqueous solution affects a wealth of structural 
and dynamic properties in a wide range of physicochem¬ 
ical, colloidal, soft-matter and biophysical systems^ . 
The standard textbook description of the electrical dou¬ 
ble layers is based on the mean-field Poisson-Boltzmann 
(PB) theory. At large surface-charge density, high 
counter-ion valency and high ion concentration - the 
so-called strong coupling limit - it is well recognized 
that PB theory fails to capture a number of qualita¬ 
tive effects, such as like-charge attraction^— and charge 
inversion^—. Liquid-state theorie s 13 ! 14 and other strong¬ 
coupling theories^^ have been employed to account for 
the strong ion-ion correlations in this regime. 

In the opposite limit - the weak-coupling regime - it 
is generally accepted that the electric double layer is well 
described by the PB theory^—. Performing a loop- 
wise perturbation expansion— in the coupling parameter 
(to be defined below), Neta^S demonstrated that the PB 
theory is the leading-order theory in the weak-coupling 
limit, and becomes exact in the limit of zero coupling 
strength. Applying Netz’s approach explicitly to surfaces 
with dielectric discontinuity, Kanduc and Podgornik 21 
concluded that, under the weak-coupling condition, the 
image force only enters as a small correction to the lead¬ 
ing PB theory, which vanishes in the limit of zero cou¬ 
pling. In particular, the self-energy due to image charge 
interaction was shown not to appear in the Boltzmann 
factor for the ion distributions. Although these demon¬ 
strations were performed explicitly for counterion-only 
systems, the conclusions are generally believed to hold 
when salt ions are added— Thus, many researchers 
in the electrolyte community consider the weak-coupling 
theory to mean the PB theory; in other words, weak cou¬ 
pling is considered synonymous with the validity of the 
PB theory. 

Physically, however, a single ion in solution next to a 
surface of a lower dielectric plate obviously should feel the 
image charge repulsion even in the absence of any sur¬ 
face charge, and the ion distribution - the probability of 
finding the ion at any location - should reflect the image 


charge interaction through the Boltzmann factor. This 
was the case studied in the pioneering work of Wagner—, 
and Onsager and Samaras^ (WOS) for the surface ten¬ 
sion of electrolyte solutions. It is rather odd that this 
interaction should become absent from the Boltzmann 
factor for the distribution of mobile ions in the weak- 
coupling limit when the surface becomes charged. It is 
also rather curious that the image interaction, which is 
absent from the Boltzmann factor in the Netz-Kanduc- 
Podgornik (NKP) approac h 15 ! 20 ' 21 in the weak coupling 
limit, “re-emerges” in the Boltzmann factor in the strong¬ 
coupling limit, though in a different form (through a fu- 
gacity expansion) 8 ' 21 ! 24 . Taking zero-surface charge as 
the limiting case of the physical weak-coupling condition, 
it is clear that the NKP and WOS approaches give dras¬ 
tically different descriptions of the same system. It is 
also difficult to physically reconcile the absence of the 
image interaction from the Boltzmann factor in the weak- 
coupling limit with its “re-emergence” in the strong cou¬ 
pling limit in the NKP approach. Furthermore, ion deple¬ 
tion near a weakly-charged dielectric interface has been 
observed in Monte Carlo simulatio n 15 ! 27 as well as pre¬ 
dicted by the hypernetted chain approximation (HNC) 
integral equation theory that includes the image charge 
interactions— 

In this work, we clarify the origin of these discrepan¬ 
cies by a re-examination of the role of the image charge 
interaction in the physical weak-coupling limit. We show 
that in the presence of a dielectric discontinuity, the phys¬ 
ical weak-coupling limit is not described by the so-called 
weak-coupling theory if the latter is meant to be the PB 
or PB with small fluctuations corrections. The image 
charge repulsion creates a boundary layer which cannot 
be captured by the the NKP approach. A nonperturba- 
tive approach yields a modified Poisson-Boltzmann equa¬ 
tion, where a screened, self-consistently determined im¬ 
age charge interaction appears in the Boltzmann factor 
for the ion concentration for any surface charge density. 
The WOS theory is an approximation of the more gen¬ 
eral framework presented here in the special case of zero 
surface charge. 

To see the origin of the boundary layer, we start by an 
analysis of the relevant length scales for the counterion- 
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only system. Consider a charged planar surface at z = 0 
with charge density a separating an aqueous solution 
(z > 0) from an semi-infinite plate (z < 0). The 
solvent and plate are taken to be dielectric continuum 
with dielectric constant eg and ep, respectively, with 
ep « es- Now consider a counterion of valency q at dis¬ 
tance 2 away from the surface. The attraction between 
the test ion and the charged surface is E sur = 2nql B az = 
z/Igg, whereas the repulsion due to its image charge is 
E im = fq 2 h/{^z), where l B = e 2 /(47 T£ 0 e s kT) is the 
Bjerrum length with e 0 denoting the vacuum permittiv¬ 
ity, Igc = l/(2irq(rl B ) is the Gouy-Chapman length and 
/ = (eg — £p)/(eg + £p) represents the dielectric contrast 
between the two media. Balancing E sur with Ei m results 
in a characteristic length: 

d=(f/2) 1/2 q(l B l GC ) 1/2 (1.1) 

Introducing the coupling parameter S = q 2 l B /lGCr — we 
see 

d ~ l B ’ET 1/2 and d/lcc ~ S 1 / 2 (1.2) 

Thus, as the coupling strength H goes to zero, d itself 
diverges, but the ratio of d to l G c (noting that Igc is 
the characteristic length scale for the double layer in the 
PB theory) goes to zero. This is a typical feature of a 
boundary layer. Physically, the competition between the 
surface charge attraction and the image charge repul¬ 
sion gives rise to a depletion boundary layer. Since the 
perturbation approach performs an expansion in powers 
of 5 15 i 20 ] 21 (which results from nondimensionalizing all 
the lengths by the longest length scale Igc), information 
within the smaller length-scale - the depletion bound¬ 
ary layer - is lost. Although this analysis is performed 
explicitly for the counterion-only system, the depletion 
boundary layer persists when salt ions are introduced. 


II. A GAUSSIAN VARIATIONAL APPROACH 

The presence of a boundary layer necessitates a non- 
perturbative treatment. Using the renormalized Gaus¬ 
sian variational approach^, one of us^H derived a general 
theory for electrolyte solutions with dielectric inhomo¬ 
geneity. In this section, we first recapitulate the key steps 
in the derivation of the general theory and then specify to 
the case of a charged plate with dielectric discontinuity. 


A. General Theory 

We consider a general system with a fixed charge dis¬ 
tribution ep ex { r) in the presence of small mobile cations 
of valency and anions of valency in a dielectric 
medium of a spatially varying dielectric function e(r). e is 
the elementary charge. The charge on the ion is assumed 
to have a finite spread given by a short-range distribution 
function h ±(r — r,) for the ith ion, with the point-charge 


model corresponding to h±( r —r^) = q±6(r — r^). The in¬ 
troduction of a finite charge distribution on the ion avoids 
the divergence of the short-range component of the self 
energy - the local solvation energy - resulting from the 
point-charge model, and reproduces the Born solvation 
energy^). However, as the emphasis of this work is on 
the long-range component of the self energy - the image 
charge interaction - which is finite for point charges, we 
will eventually take the point-charge limit for the ion. 
The diverging but constant local solvation energy in the 
point-charge limit can be regularized by subtracting the 
same-point Green function in the bulk, as discussed be¬ 
low. Since we work in the low concentration regime for 
the ions (c < 0.1 M) (the Debye-Hiickel regime), the ex¬ 
cluded volume effects of the ions are unimportant, and 
so we treat the ions as volumeless. 

The total charge density including both the fixed 
charges and mobile ions is 

ep{ r) = e p ex { r) + e J dr'[h+(r' - r)c + (r') 

— h-(r' — r)c-(r')] (2.1) 

with c±(r) = Y^i=i H r ~ r i) the particle density oper¬ 
ator for the ions. The Coulomb energy of the system, 
including the self energy , is 

H = e — j drdr'p(r)G 0 (r,r')p(r') (2.2) 

where Go(r, r') is the Coulomb operator given by 

- V • [e 0 e(r)VG 0 (r, r')] = 5(r - r') (2.3) 


It is convenient to work with the grand canonical parti¬ 
tion function 
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where p± are the chemical potential for the cations and 
anions, and v± are some characteristic volume scales, 
which have no thermodynamic consequence. We perform 
the usual Hubbard-Stratonovich transformation to Eq. 
12.41 by introducing a field variable </>(r), which yields 


fl = 


Zn 


Dcfex p {— L[<fi]} 


(2.5) 


The “action” L is 

L = j dr [ le(V^) 2 + ip ex (j) - TA+e"^ 

-TA_e^-^] (2.6) 

Z 0 is the normalization factor given by 


^o = 


Dip exp 



dre(V</>) 2 


[det G 0 ] 1/2 (2.7) 
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where Gq 1 = V r • [e(r)V r /<5(r — r')] is the inverse of the 
Coulomb operator in Eq. 12.31 e = e 0 e/(^e 2 ) is a scaled 
permittivity, and \± = e ll± /v± is the fugacity of the ions. 
We have used the short-hand notation h±<j> to represent 
the local spatial averaging of <p by the charge distribution 
function: h±<j) = f dr'h±(r' — r)tp(r'). The function T(r) 
in Eq. 12.61 is introduced to constrain the mobile ions to 
the solvent region. 

Equations 12. 5l and 12.61 are the exact field-theoretic rep¬ 
resentation for the partition function. Because the action 
is nonlinear, the partition function cannot be evaluated 
exactly. The lowest-order approximation corresponds to 
taking the saddle-point contribution of the functional in¬ 
tegral, which results in the Poisson-Boltzmann equation. 
A systematic loop expansion can be developed to account 
for fluctuations around the saddle point in an order by 
order manner. In practice, most theoretical treatments 
only include one-loop corrections. The loop expansion 
involves expanding the action around the saddle point in 
polynomial forms. However, the fluctuation part of the 
electrostatic potential due to image charge interaction 
becomes very large near the dielectric interface; thus any 
finite-order expansion of the (e Tlh±< ^) term, which be¬ 
comes the Boltzmann factor in the ion distribution (see 
Eq. 12.111) . is problematic. The absence of the image- 
charge self-energy in the Boltzmann factor in the pertur¬ 
bation approachesi^J&^lT— is thus a consequence of low- 
order expansion of the exponential function of an imagi¬ 
nary variable when the variable can be quite large. 

To develop a nonperturbative theory, we perform a 
variational calculation of Eq. 12.51 using the Gibbs- 
Feynman-Bogoliubov bound for the grand free energy 
W = — lnfi, which yields 

W<-lnn ref + (L[<p\-L ref [<Pl) (2.8) 

where 

Mref = vr [ D(j )exp {-L ref [0]} (2.9) 

The average (• • •) is taken in the reference ensemble with 
the action L re f. We take the reference action to be of 
the Gaussian form centered around the mean —iip 

Lref = \ j drdr'[(/)(r) + ^(r)]G _1 (r, r')[<p{r') + iip{ r')] 

( 2 . 10 ) 

where G 1 is the functional inverse of the Green func¬ 
tion G, and the introduction of i is to keep the mean 
electrostatic potential ip real, ip and G are taken to be 
variational parameters for the grand free energy func¬ 
tional. 

With the Gaussian reference action Eq. 12.101 all the 
terms on the r.h.s. of Eq. 12.81 can be evaluated analyti¬ 
cally (see Appendix A for detailed derivation). The lower 
bound of the free energy is obtained by extremizing the 
r.h.s. of Eq. 12.81 with respect to ip and G, which results 


in the following two variational conditions: 

- V • (eViP) = p ex + T\ + q + e-^- u + - rA_t7_e^-“- 

(2.11) 

- V ■ [eVG(r, r')] + 2/(r)G(r, r') = <5(r - r') (2.12) 

where u± is the self energy of the ions 

u±(r) = - J dr'dr"h±(r — r')G(r',r")h±(r" — r) (2.13) 

/(r) = [g 2 c+(r) + qf c_(r)] /2 is the local ionic strength, 
with the concentration of cations and anions given by 

c±(r) = A±Texp [=R±lH r ) ~ u ±( r )} ( 2 - 14 ) 

Eqs. 12.11112. In forms a set of self-consistent equations 
for the mean electrostatic potential ip( r), the correlation 
function (Green function) G(r, r') and the self energy 
u± (r) of the ions, which are the key equations for weakly 
coupled electrolyte s 30 ’ 32 . Eq. I2.11l has the same form as 
the PB equation, but now with the self-energy of the ions 
appearing in the Boltzmann factor. The appearance of 
the self energy in the Boltzmann factor reflects the non¬ 
linear feedback of the fluctuation effects, an aspect that 
was missing in a perturbation expansion. The self-energy 
given by Eq. I2.13l is a unified expression that includes the 
Born energy of the ion, the interaction between the ion 
and its ionic atmosphere, as well as the distortion of the 
electric field by a spatially varying dielectric function, the 
latter taking the form of image charge interaction near 
a dielectric discontinuity. In general, the self energy is 
spatially varying if there is spatial inhomogeneity in ei¬ 
ther the dielectric constant or the ionic strength. Making 
use of the variational conditions Eqs. I2.11l and l2.12l and 
evaluating the fluctuation part of the free energy arising 
from Gaussian integrals by using the charging method (as 
shown in Appendix B), we obtain a simple expression for 
the equilibrium grand free energy: 

W = - J dr (c+ + c_) + i J drip ( p ex - q+c+ + q-cJ) 
+ J drl(r) J dp [G(r, r; p) — G(r, r)] (2-15) 

where p is a “charging” variable. G(r, r; p) is the same- 
point Green function obtained from solving Eq. I2.12l but 
with the term /(r) replaced with pl( r). Note that the 
free energy expression Eq. I2.15l is finite even in the point- 
charge limit. Although both G(r, r; rj) and G(r, r) are in¬ 
finite, their divergent parts exactly cancel; the remaining 
difference is finite and accounts for the leading-order ion- 
ion correlation effect. Unlike previous field-theoretical 
treatment s 20 : 22 : 31 , no microscopic cut-off is needed in our 
theory. 

B. Weakly Charged Plate 

We now specify to the case of a charged plate with 
dielectric discontinuity in contact with an electrolyte so¬ 
lution. The fixed external charge density is then p ex (r) = 
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aS(z). For concreteness, we take the surface charge to be 
positive. Both T and e(r) are now step functions: T = 0 
and e(r) = ep for 2 < 0; F = 1 and e(r) = eg for 2 > 0. 
In the solvent region (z > 0), Eq. 12. Ill becomes 

- = \+q+e-^- u + - A _q_e^- u - (2.16) 


with the boundary condition {dij)/dz) z - o = — a/es , 
which is obtained by integrating Eq. I2.11l between 2 = 0“ 
and z = 0 + and noting that {d^/dz) = 0 for 2 < 0. Since 
the solvent has a uniform dielectric constant, the Born 
energy is constant and can be absorbed into the refer¬ 
ence chemical potential. It is then convenient to single 
out this constant contribution by rewriting Eq. 12.131 as 

u±{v) = i J dv'dv"h ± (v - r')^— A— - r) 

+ — [ dv'dv"h±(v — r') G( r',r")-:-- 

2 J ±K 47re s |r'-r"|_ 

xh ± { r"-r) (2.17) 

The first term gives a constant Born energy of the ion 
q±/8iresa±, with a± the Born radius of the ion^S. The 
remaining term is finite in the point-charge limit. We 
can thus take the a± —> 0 limit for this term in the fi¬ 
nal expression, or equivalently and more conveniently by 
directly take the point-charge limit in the distribution 
h±( r — r') = <5(r — r'). The nontrivial and nonsingular 
part of the self energy u± is then 


* 4 ( r ) = % 1 ™ 

Z i*'— 


G(r,r') - 


47res|r — r'| 


(2.18) 


To avoid the complexity of solving the equation for the 
Green function (12.1211 , previous work usually invoked ap¬ 
proximate schemes, e.g., by replacing the spatially vary¬ 
ing screening length by the bulk Debye lengt h 14 ’ 26 ’ 27 ’ 35 ’ 36 
or using a WKB-like approximation^— . However, the 
screening on the image force at the dielectric interface 
is inhomogeneous, long-ranged and accumulative, which 
cannot be captured fully by these approximate methods. 
In this work, we perform the full numerical solution of 
the Green function, which provides the most accurate 
treatment of the inhomogeneous screening effect at the 
dielectric interface. To solve the Green function in the 
planar geometry, it is convenient to work in a cylindrical 
coordinate system (r, z). Noting the isotropy and trans¬ 
lational invariance in the directions parallel to the sur¬ 
face, we can perform a Fourier transform in the parallel 
direction to write: 

1 r°° . 

G(r,z,z')=— / kdkJ 0 (kr)G(k, 2 , z’) (2.19) 

2tt Jo 

where Jo is the zero-order Bessel function. G(fc, 2 , z') 
now satisfies: 

- d ^ + [k 2 (z) + k 2 ] G(k, 2 , z') = —(5(2, z’) 

oz z es 

( 2 . 20 ) 


for 2 > 0, with the boundary condition esdG/dz — 
ktpG = 0 at 2 = 0. k(z) = [2I(z)/es] 1 ^ 2 can be con¬ 
sidered the inverse of the local Debye screening length. 

Eq. 12.201 is solved numerically by using the finite differ¬ 
ence metho d 37 ’ 38 . The free-space Green function satisfy¬ 
ing —d 2 Go/dz 2 + k 2 Go = d(z,z')/es, though analytically 
solvable, is also solved numerically along with Eq. 12.201 
to ensure consistent numerical accuracy in removing the 
singularity of the same-point Green function. The non- 
divergent part of the self energy is then: 


*4(2) 



G(k, 2 , z) - G 0 (k, 2 , z) 


kdk 


( 2 . 21 ) 


Far away from the plate surface (2 —> 00 ), the ion con¬ 
centration approaches the bulk value c^_, and from Eq. 
12.141 (where we set ipb = 0, or equivalently absorbing a 
constant </>& into the definition of the fugacity), the fugac- 
ity of the ions is given by X± = c± exp [— q±Kb/(Sires)] 
where Kb is the inverse screening length in the bulk. Note 
that this relationship automatically takes into account 
the Debye-Huckel correction to fugacity due to ion-ion 
correlations. 

The theory presented above is derived explicitly with 
added salt. However, application to the counterion- 
only system is straightforward through an ensemble 
transformation^. 


III. NUMERICAL RESULTS AND 
DISCUSSIONS 

In this section, we apply the theory presented in the 
last section to an electrolyte solution in contact with a 
weakly charged plate. We first examine the counterion- 
only system to highlight the depletion boundary layer 
issue and then study the consequences of the deple¬ 
tion boundary layer on the structure and thermodynamic 
properties of the electric double layer with added salts. 


A. Counterion-only Case 

For the counterion-only system, the PB theory ad¬ 
mits an analytical solution for the counterion distribu¬ 
tion: c(z) = 1/ [2nlBq 2 {z 2 + Iq C )\, which is character¬ 
ized by a single length scale, the Gouy-Chapman length. 
The counterion concentration profile is shown in Fig. 1 as 
the dashed line, which decays monotonically. In contrast, 
when there is dielectric discontinuity, our theory predicts 
a qualitatively different behavior. The presence of the de¬ 
pletion boundary layer inside the Gouy-Chapman length 
is obvious, and is consistent with results from Monte 
Carlo simulatio n 15 ’ 27 . Within the depletion boundary 
layer (2 < d ~ ZbS - 1 / 2 ), image charge repulsion is domi¬ 
nant and ions are excluded from the plate surface. In the 
point-charge model, the self energy diverges to infinity 
at the plate surface; thus the ion concentration vanishes 
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at z = 0. The vanishing of the ion concentration obvi¬ 
ously contracts the PB prediction but is also incapable of 
being captured by any perturbation corrections around 
the PB limi t 1 7 i 18 i 20 ~ —. Simply put, these perturbative 
approaches fail to satisfy the boundary condition for the 
ion concentration at z = 0, as is typical with boundary 
layer problems. Beyond the depletion boundary layer 
(z > d), surface charge attraction prevails and the ion 
concentration approaches the PB profile sufficiently far 
away from the surface. 



z[nm| 

FIG. 1. Ion concentration for the counterion-only system 
for different ion valency q. es = 80 and ep = 2.5. a = 
l/(100q)(e/nm 2 ). The Gouy-Chapman length is kept con¬ 
stant ( Igc = 22.7 nm) for counterions of different valencies. 
The coupling parameter E is 0.031q 2 . 

The PB theory predicts a universal profile q 2 c(z) 
for counterions of different valencies when the Gouy- 
Chapman length is kept fixed. However, from our scal¬ 
ing analysis in the Introduction, the depletion boundary 
layer should increase linearly with valency; see Eq. 11.11 
This prediction is borne out by our numerical result as 
shown in Fig. 1. Therefore, the boundary layer problem 
becomes more severe for ions of high valency. The scaling 
of the boundary thickness with the coupling parameter 
predicted from Eq. I1.2l is also confirmed by our numerical 
results (data not shown). 


B. With Symmetric Salt 

When there are added salt ions in the solution, the 
image force affects the distribution of both the counte¬ 
rions and coions. The PB theory predicts that the dou¬ 
ble layer structure is characterized by the Debye screen¬ 
ing length k - 1 under the condition that kT 1 -C Igc , 
with a monotonically decreasing counterion and mono- 
tonically increasing coion distribution. In contrast, both 
the counterion and coion concentration must vanish at 
the surface, but their approach to the bulk concentration 
is different: the coions increases monotonically, while the 


counterions goes through an overshoot. Furthermore, we 
find two regimes depending on the relative width of the 
screening length and the boundary layer thickness, which 
itself is in turn affected by the screening. At low salt con¬ 
centration, kT 1 d and ion depletion is confined in a 
boundary layer very close to the plate surface; both the 
ion distribution and electrostatic potential approach the 
profile predicted by PB beyond the boundary layer. As 
the salt concentration increases, the width of the deple¬ 
tion boundary layer becomes comparable to the screen¬ 
ing length and the two length scales remain comparable 
thereafter; the image charge interaction then affect the 
entire range of the double layer. In Figure 2 we show 
the ion distribution of a 0.1M 1:1 electrolyte calculated 
by our theory. The contrast with the PB result is quite 
striking. 



FIG. 2. Ion concentration profile for a 1:1 electrolyte solution 
with c b = 0.1 M. es = 80, ep = 2.5 and a = le/100 nm 2 . 

The change in the double layer structure will affect 
a wealth of interfacial properties. As an example, we 



FIG. 3. The surface energy f s as a function of the salt con¬ 
centration for a 1:1 electrolyte solution, es = 80, ep = 2.5 
and a = le/100 nm 2 . 
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show in Figure 3 the surface excess free energy f s = 
f 0 (w — w b )dz (where w is the grand free energy density 
and w b is its bulk value) as a function of the salt concen¬ 
tration. The PB theory predicts a monotonic decrease of 
f s that scales approximately with (c b ) -1 / 2 , which arises 
from the electric field contribution in the free energy due 
to the surface charg o 39 ! 40 . With the inclusion of image 
charge interaction, our theory shows that f s changes non- 
monotonically. At low salt concentration (c b < 10~ 3 M), 
f s calculated by our theory follows closely the PB result; 
this is because the region affected by the image charge 
repulsion is relatively narrow compared to the screening 
length, giving a relatively small contribution to the sur¬ 
face excess energy when integrated over the entire solu¬ 
tion. As the salt concentration increases (c b > 10~ 2 M), 
our theory predicts a sign change in the slope of f s vs. 
c b : f s increases with increasing c b , opposite to the PB re¬ 
sult. In this concentration regime, the width of the deple¬ 
tion boundary layer is comparable to the Debye screening 
length, and the entire double layer region is affected by 
the image charge interaction as shown in Figure 2. The 
increase in f s is now largely due to the depletion (i.e., 
negative adsorption) of mobile ions. The slope of log(/ s ) 
vs log(c b ) is less than 1 because of the increased screen¬ 
ing of the image force as the salt concentration increases. 
The sign change of df s /dc b corresponds to the crossover 
in the length scale relationship from k -1 3> d to re -1 ~ d. 
As the excess surface energy determines the spreading of 
a liquid drop on a solid surface, this result implies a qual¬ 
itatively different behavior for the spreading of a drop of 
electrolyte solution than that predicted by the PB theory. 
We also note that the nonmonotonic behavior discussed 
here shares the same physics as the Jones-Ray effect^— 
for the interfacial tension observed at the water/air and 
water/oil interfaces. 


C. With Asymmetric Salt 

The effects of image charge become more complex if the 
salt ions are of unequal valency. Because of the quadratic 
dependence of the image force on the valency, the higher- 
valent ions are pushed further away from the surface, ne¬ 
cessitating a compensation by the lower-valent ions in the 
space in between. The difference in the image force be¬ 
tween the counterions and the coions induces additional 
charge separation and hence electric Held within the de¬ 
pletion boundary layer. The induced net charge within 
the boundary layer alters the effective surface charge, 
which can affect the double layer structure outside the 
boundary layer. For the case where the coions are of 
higher valency than the counterions, the induced electric 
field due to unequal ion depletion counteracts the field 
generated by the surface charge. With the increase of 
the salt concentration, the induced field can exceed that 
generated by the bare surface charge, leading to a sign 
change in the effective surface charge known as charge 
inversion. The double layer structure becomes qualita¬ 


tively different from that predicted by the PB theory as 
shown in Figure 4: the electrostatic potential is of the 
opposite sign to the PB result. Excess counterions ac¬ 
cumulate in the depletion boundary layer, overcharging 
the plate surface, while the coions are enriched outside 
the boundary layer, serving to screen the inverted surface 
charge. In this case, the PB theory qualitatively fails to 
describe the entire double layer structure. 



FIG. 4. Charge inversion for a 0.05M 2:1 electrolyte solution 
near a positively charged plate, (a) Dimensionless electro¬ 
static potential and (b) net charge density (<j+c+ — g_c_). 
Es = 80, Ep = 2.5 and a = le/100 nm 2 . 


D. Uncharged Surface: Image Charge vs. 

Correlation Effect 

The case of an electrolyte solution next to an un¬ 
charged surface (cr = 0) reduces to the problem treated 
by Wagner, Onsager and Samaras. The self energy due 
to image charge repulsion appears in the Boltzmann fac¬ 
tor and is responsible for the depletion layer in the ion 
distribution near the surface as shown in Figure 5. Note, 
however, in the original WOS theory as well as in subse¬ 
quent treatment s 14 ! 26 ' 27 ! 35 ’ 36 , the image charge term was 
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added to the Boltzmann factor ad hoc based on physical 
intuition, whereas in our theory, its appearance is the re¬ 
sult of systematic derivation. Therefore, our theory not 
only recovers the WOS theory (upon making additional 
approximations, e.g., by using the constant bulk screen¬ 
ing length for the image force potential) but also provides 
the means for systematically improving the WOS the¬ 
ory. First, our theory captures the anisotropic screening 
cloud around an ion near the interface due to the spatially 
varying ion concentration near the surface. The inhomo¬ 
geneous ionic cloud in the depletion layer and its effect on 
the screening of the test ion are treated self-consistently 
in our theory, whereas this inhomogeneous screening is 
missing in the WOS theory. Second, by including the 
mean electrostatic potential generated by charge separa¬ 
tion, our theory can describe salt solutions with unequal 
valency such as the case of 2:1 electrolyte shown in Figure 
5b. Finally, our theory provides a more accurate expres¬ 
sion for the excess free energy by properly accounting for 
the inhomogeneous screening effect and the fluctuation 
contribution to the free energy. Thus, we expect our the¬ 
ory to be able to better predict the surface tension of 
electrolyte solutions in comparison to the WOS theory, 
especially at higher salt concentrations (where accurate 
treatment of the screening becomes more important). 

The inhomogeneous screening results in a correlation 
effect that can lead to ion depletion near the surface^: an 
ion interacts more favorably with its full ionic atmosphere 
far away from the surface than in the vicinity of the sur¬ 
face. This correlation effect is stronger for multivalent 
ions, which pushes them further away from the interface 
than the monovalent ions. The correlation-induced ion 
depletion near the surface can take place both with and 
without the dielectric contrast, and is well captured by 
our theory, as shown in Figure 5. While the ion depletion 
without the dielectric contrast is induced by the corre¬ 
lation alone, the ion depletion in the presence of the di¬ 
electric contrast is due to both the correlation effect and 
the image charge effect, which enhance each other. As a 
result, both ion depletion, as well as charge separation in 
the case of 2:1 electrolyte, are more pronounced in the 
presence of image charge than due to correlation alone. 

Ion depletion due to correlation alone is most notice¬ 
able when the surface is uncharged. When the surface 
is charged, the surface attraction for the counterions will 
dominate over such correlation effect in the absence of 
image charge repulsion. In contrast, depletion due to 
image charge repulsion persists for both the counterions 
and coions even when the surface is charged. 


IV. CONCLUSIONS 

In this work, we have shown that the image charge 
repulsion creates a depletion boundary layer near a di¬ 
electric surface, which cannot be captured by a regular 
perturbation method. Using a nonperturbative approach 
based on a variational functional formulation, we find 




FIG. 5. (Color online) Ion concentration (scaled by the bulk 
ion concentration c±) for (a) 0.01M 1:1 electrolyte solution 
(c+ = C-) and (b) 0.01M 2:1 electrolyte solution near an 
uncharged interface (o = 0) with dielectric contrast (es = 
80, ep = 2.5) in comparison with the case without dielectric 
contrast (es = ep = 80). Profiles calculated by our theory 
are shown by colored lines; results from the PB theory are 
given as black dot lines. 


that the self energy of the ion, which includes contribu¬ 
tions from both the image charge interaction and the in¬ 
terionic correlation, appears explicitly in the Boltzmann 
factor for the ion distribution, resulting in a self-energy 
modified Poisson-Boltzmann equation as the appropriate 
theory for describing the physical weak-coupling condi¬ 
tion. This image-charge self energy is not diminished by 
reducing the surface or the ionic strength in the solution; 
in the presence of a significant dielectric discontinuity, 
there is no limiting condition for which the PB theory 
is valid. For zero surface charge, our theory reduces to 
the WOS theory upon further approximations. Thus, our 
theory provides both the justification for the WOS theory 
and means for systematically improving the WOS theory, 
for example, by including the mean electrostatic potential 
















generated by the charge separation in salt solutions with 
unequal valency or other asymmetries between the cation 
and anions, such as different size and polarizability^. 

The weak-coupling condition in the presence of dielec¬ 
tric discontinuity covers many soft-matter and biophysi¬ 
cal systems. Many phenomena, such as the surface ten¬ 
sion of electrolyte solution s 43 ! 44 , salt effects on bubble 
coalescence^, and the ion conductivity in artificial and 
biological ion-channels^—, cannot be explained, even 
qualitatively, by the PB theory. The presence of the im¬ 
age charge interaction results in a very different picture 
of the electrical double layer from that provided by the 
PB theory, and can give rise to such phenomena as like- 
charge attraction and charge inversion even in the weak- 
coupling conditions*; these phenomena have usually been 
associated with the strong-coupling condition. The PB 
theory has played a foundational role in colloidal and in¬ 
terfacial sciences: the DLVO theory, interpretation of the 
zeta potential, experimental determination of the surface 
charge and the Hamaker constant, are all based on the 
PB theory^. With the inclusion of the image charge in¬ 
teraction, some of the well known and accepted results 
will have to be reexamined. 


and 


( e =F<&±x) = exp [ - i J dr'dr"h±(r - r')G(r',r") 

x h±(r" — r)] (A4) 

Substituting Eqs. IA21IA4I into Eq. IA1I we obtain the 
variational form of the grand free energy as 


W= 2 h 1 


det G 
det G 0 


~\J rfre(V ^) 2 


~\J dvdv' [G _1 (r, r') — Gg 1 (r, r')] G(r,r') 

+ J dr ( p ex il> - PA + e-^-“+ - TA_< 




)(A5) 


where u± is the self energy of the ions given by Eq. 12.131 
Minimizing Eq. lASI with respect to ip and G gives rise to 
Eq. I2.11l and Eq. I2.12l in the main text. 


Appendix B: Simplification of the fluctuation 
contribution in the free energy 
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Appendix A: Derivation of the key equations in 
Section II.A 

We define x = <P + iip as the fluctuation part of the 
field (p, which is a Gaussian variable by our ansatz. The 
variational grand free energy can be approximated by the 
r.h.s of Eq. [Has 


W = — J dr [c+ (r) + c_ (r)] 

+ \ j drip(r) [p ex ( r) - q + c+( r) + g_c_(r)] 

-^"(s^)-/* /(r)G,(r - r > (B1) 

The last two terms in Eq. IB II are due to the fluctuation 
contribution. We note that the Green function equation 
(Eq. 12.121) can be written in the matrix form as 

Go 1 G + 2/(r)G =1 (B2) 


W — lnf2 re / + (L [cp] L re f [(/>]} — ^ l n q 0 ^) 

- \ J drdr'{6(r' - r) [e(VV’) 2 - e((Vy) 2 )] + G _1 (r, r') 
x <x(r)x(r , ))} + J dr[p ex iP-r\ + e- q ^{e~ ih +x) 


where I is the identity matrix (not to be confused with 
the local ionic strength I). Right multiplication of the 
above equation by G~ 4 , we obtain 

G” 1 = Go 1 + 2/(r)I (B3) 

Note also that 


+ rA_e 9 ~^(e^ -x )] (Al) 

The averages in Eq. lAll can be evaluated exactly because 
the distribution of x is Gaussian. Noting that 

(x(r)x(r')) = G(r, r') (A2) 


In 



det G A 
det Go/ 



= In det G — In det Go 


6 In det G 
(5G _1 (r, r') 


5G-\ r, 


r') 


(B4) 


we have 


drdr'Spr 


rWxf) 


The innermost integral is a functional integration over 
G _1 from G,/" 1 to G _1 . Since In det G is the result of a 
Gaussian functional integral, we have 


drdr'Vr • [e(r)V(.<5(r — r')] G(r, r') (A3) 


S In det G 
<5G -1 (r, r') 


(B5) 


G(r,r') 
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Therefore, 



(B6) 


As the integration goes from G 0 1 to G 1 , the integrand 
changes from Go to G. From Eq. lB3l a convenient 
path for integrating Eq. IBfil is to introduce a continuous 
“charging” variable 77 that goes from 0 to 1 multiplying 
the 2 1 term in Eq. IB31 while keeping the density profile 
fixed. Obviously the Green function is Go for 77 = 0 and 
is G for rj = 1. For any intermediate value, we denote the 
Green function as G(r, r';r)). Using Eq. IB31 the above 


integral becomes, 

h, {^) = - 2 J dr J dr 'l /(r) ' S(W) 

x G(r, r'; r))drj = —2 J drl(r) j G(r, r; 77)^77 (B7) 

where G(r, r; 77 ) is to be understood as the limit 
G(r, r; 77 ) = linir'-^r G(r, r'; 77 ), and the Green function 
G(r, r'; 77 ) is the solution of 


- V ■ [eVG(r, r')] + 2r)I(r)G(r, r') = S(r - r') (B 8 ) 

With Eq. IB71 the fluctuation contribution to the free 
energy is 


1 / det G 

- In 


+ J dr/(r)G(r,r) 


z \ det Go, 

= - j dr Hr) J dr) [G(r, r; 77 ) - G(r, r)] (B9) 

which is finite even in the point-charge limit. 
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